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1.  Introduction 

The  problem  we  consider  is  the  solution  of  the  linear  system 

Ax  =  6,  (1.1) 

where  A  €  is  symmetric  positive  definite  but  possibly  nearly  singular.  We  will  denote  the 

eigenvalues  of  A  by  Aj  <  A2  <  .  ..,<  A and  the  corresponding'orthonormal  eigenvectors  by 
{u>i,  u>2  -  -  • ,  Wtf}.  For  simplicity,  we  will  assume  that  the  near-nullity  of  A  is  at  most  one,  i.e., 
that  no  eigenvalue  other  than  A|  could  be  close  to  zero.  The  positive  definiteness  assumption  of  A 
can  be  slightly  relaxed,  by  replacing  it  with  the  less  restrictive  assumption  that  a  few  of  the  first 
eigenvalues  of  A  are  negative. 

The  solution  to  (1.1)  can  be  expressed  in  the  form: 

wTb 

x  =  xd  +  -A—  t»i,  (1.2) 


where 


_  A  xvjb 

~t~w' 


In  the  above  expression,  we  have  separated  the  part  of  x  in  the  direction  of  wi  from  the  part  Xj 
that  is  orthogonal  to  it.  The  vector  X4  is  called  the  deflated  solution  of  (1.1)  and  (1.2)  is  called 
the  deflated  decomposition  of  x.  When  At  is  small  but  wfb  is  not  small,  then  the  last  term  in 
(1.2)  will  dominate  xj,  and  in  finite  precision  arithmetic,  it  would  be  difficult  to  recover  x&  from 
x  with  close  to  full  machine  precision.  Therefore,  when  A  is  nearly  singular,  it  is  often  more 
appropriate  to  compute  the  deflated  decomposition  of  x  rather  than  to  compute  x  directly.  This 
requires  computing  Xd,  Ai  and  u>\. 

The  deflated  decomposition  is  useful  in  many  applications,  such  as  when  solving  bordered 
singular  systems  [5,  4,  1]  that  arise  in  continuation  methods  for  solving  nonlinear  systems  [3,  9] 
and  bifurcation  computations  [8,  15]  and  in  constrained  optimization  problems[7]. 

Of  course,  one  could  compute  the  deflated  decomposition  by  first  computing  the  eigenvalue 
decomposition  of  A,  but  this  is  often  too  expensive  for  large  problems.  In  earlier  work  by  Stewart 
[20]  and  Chan  [2],  the  deflated  decomposition  is  computed  by  implicit  algorithms  in  which  only  a 
direct  solver  for  A  (such  as  an  LU  factorization  method)  is  needed.  In  this  paper,  we  look  at  how 
iterative  methods  based  on  the  Lanczos  process  can  be  used  to  compute  the  deflated  decomposition 
in  a  numerically  stable  manner.  The  objective  is  to  be  able  to  compute  the  deflated  decomposition 
by  only  accessing  A  in  the  form  of  a  matrix-vector  product  Av.  Such  a  method  has  obvious 
advantages  for  large  problems.  Basically,  the  Lanczos  method  produces  a  small  tridiagonal  matrix 
T  approximately  similar  to  A  and  the  deflation  techniques  in  [2,  20]  can  then  be  applied  to  T.  In 
addition,  we  propose  a  new  method  based  on  the  QL  iteration  for  deflating  the  tridiagonal  matrice 
T. 

Note  that,  as  explained  above,  the  naive  method  of  applying  a  Lanczos  type  procedure,  such 
as  the  symmetric  conjugate  gradient  method,  to  (1.1)  directly  will  fail  to  compute  14  accurately 
when  A  is  nearly  singular.  If  w\  is  known  a  priori ,  then  a  modified  conjugate  gradient  algorithm, 
in  which  the  iterate  at  each  iteration  is  orthogonalized  with  respect  to  [10],  might  be  expected 
to  be  effective.  The  algorithms  to  be  presented  in  this  paper  do  not  require  a  priori  knowledge  of 
w j  or  A|. 

In  Section  2,  we  review  the  deflation  techniques  of  [2,  20]  and  present  the  new  QL  algorithm 
for  tridiagonal  matrices.  In  Section  3,  we  review  the  Lanczos  process  and  then  in  Section  4,  we 
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explain  how  to  combine  the  techniques  of  Sections  2  and  3  to  compute  the  deflated  decomposition 
of  (1.1).  In  Section  5,  we  discuss  some  implementation  issues  concerning  the  loss  of  orthogonality 
of  the  Lanczos  vectors  and  stopping  criteria.  In  Section  6,  we  present  a  convergence  analysis  that 
essentially  shows  that  the  deflated  solution  converges  at  the  rate  that  would  have  been  achieved  if 
we  removed  the  eigenvalue  Ai  from  the  spectrum  of  the  original  matrix  A.  Finally,  we  present  some 
numerical  results  in  Section  7  and  end  in  Section  6  with  some  concluding  remarks.  Throughout 
this  paper,  upper  case  Latin  letters  denote  matrices,  lower  case  Latin  letters  denote  vectors  and 
lower  case  Greek  letters  denote  scalars.  We  will  use  the  notation  ||  •  ||  to  denote  the  2-norm. 

2.  Deflated  Decomposition  for  Tridiagonal  Systems 

In  this  section,  we  discuss  techniques  for  computing  the  deflated  decomposition  of  solutions  to 
the  linear  system 

Tz  =  /,  (2.1) 

where  T  is  tridiagonal.  First,  we  review  the  deflation  techniques  of  [2,  20],  which  are  based  on 
orthogonal  projections.  These  techniques  are  designed  for  general  linear  systems  but  can  be  applied 
to  tridiagonal  systems  to  produce  an  efficient  deflation  algorithm.  Next  we  will  present  a  new 
deflation  algorithm  based  on  the  QL  iteration  specifically  designed  for  tridiagonal  matrices. 

2.1.  Deflation  by  Orthogonal  Projection 

Let  P  =  /  —  u\ u[  denote  the  orthogonal  projector  with  respect  to  the  eigenvector  ui  cor¬ 
responding  to  the  smallest  eigenvalue  Ai  of  T.  Then  the  deflated  solution  Zd  of  Tz  =  f  can  be 
characterized  as  the  unique  solution  to  the  following  singular  but  consistent  system  with  a  con¬ 
straint: 

PTzd  =  Pf 
Pzd  =  Zd • 

Based  on  this  characterization,  Stewart  [20]  and  Chan  [2]  propose  the  following  implicit  algorithm 
for  computing  Zd . 


Algorithm  Deflate: 

1.  Compute  Ai  and  «i  of  T  by  a  few  steps  of  inverse  iteration. 

2.  Solve  the  system  Tz  =  Pf  for  z. 

3.  Compute  Zd  =  Pz. 

It  is  shown  in  [2]  that  Algorithm  Deflate  computes  Zd  in  a  stable  manner.  For  tridiagonal 
matrices,  the  above  algorithm  can  take  full  advantage  of  efficient  tridiagonal  solvers. 

2.2.  Deflation  by  the  QL  Method 

Assume  that  we  have  computed  the  eigenpair  (Ai,  uj)  of  the  tridiagonal  matrix  T  by  inverse 
iteration.  The  main  idea  behind  the  QL-deflation  method  is  that  if  we  apply  one  step  of  the  QL 
iteration  [14]  for  computing  the  eigenvalues  of  T  with  the  shift  Ai,  then  the  resulting  transformed 
tridiagonal  matrix  decouples  in  a  way  which  allows  the  deflated  decomposition  to  be  computed 
easily  and  in  a  stable  way. 

Specifically,  the  QL  transformation  amounts  to  first  computing  the  LQ-factorization  of  T-Aj/ , 


T- Xil=  LQ, 


(2.2) 


where  L  is  lower  triangular  and  Q  is  orthogonal,  and  performing  the  product  in  reverse  order, 
adding  back  the  shift: 


r(l>  =  QL  +  A,/. 


(2.3) 


The  transformation  from  T  to  is  the  well  known  QL-transformation  with  shift  Ai  [I4j.  In 
practical  implementations  the  two  operations  (2.2)  and  (2.3)  are  performed  in  one  single  pass.  It 
can  be  shown  that  T M  has  the  form  shown  below 


A 1 1  0  0 
0 
0 

r(‘)=  f 


where  T  is  tridiagonal.  It  is  well  known  that  is  unitarily  similar  to  T  since 

T*1)  =  QL  +  A  J  m  Q[T  -  \iI)Qt  +  A X1  m  QTQt. 

The  system  Tz  =  /  can  now  be  transformed  into  one  for  T^\  namely, 

TWy  =  Qf  =  f,  (2.4) 

where  y  —  Qz.  If  we  partition  y  and  /  into  (yi,yi)T  and  ( f\,h)T  to  conform  with  the  partition  of 
rW,  then  the  deflated  solution  yd  of  (2.4)  is  easily  seen  to  be: 

y«  =  (o,r -lh)T. 

The  deflated  solution  Z4  of  (2.1)  can  then  be  obtained  by 

**  =  QTVd- 

Since  T  is  nonsingular  by  the  assumption  that  the  nullity  of  A  is  at  most  one,  this  approach 
requires  only  solutions  of  nonsingular  tridiagonal  systems. 
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3.  The  Lanczos  Algorithm  for  Solving  Linear  Systems. 

In  this  section  we  put  aside  the  issue  of  near  singularity  temporarily  and  describe  the  Lanczos 
method  for  solving  linear  systems  [13].  Consider  the  linear  system: 

Ax  =  b  (3.1) 

Suppose  that  a  guess  x to  the  solution  is  available  and  let  tq  be  its  residual  vector:  tq  =  b  —  Ax^ . 
Then  the  Lanczos  algorithm  for  solving  (3.1)  can  be  described  as  follows: 


Algorithm:  The  Lanczos  algorithm  for  solving  linear  systems 

(1)  Stage  1  :  Generate  the  Lanczos  Vectors 

•  Start:  v\  =  r0/||ro|| 

•  For  j  =  1, 2 . . .  m  compute 

Qj  {Avj,vj)  (3.2) 

vj+i  :=  Avj  -  otjVj  -  0jVj.  i,  (0i  vo  =  0)  (3.3) 

0j+l  ■=  11*7+1  II  (3.4) 

vj+ 1  :=  Vj+i/0j+i  (3-5) 

(2)  Stage  2  :  Form  the  approximate  solution 

a.(m):=a.(°)  +  Vmr-1(||ro||ei)  (3.6) 

where  Vm  =  [vx,  i>2,  ,.t)m]  and  Tm  is  the  tridiagonal  matrix  Tridiag[0j,ctj,0j+1]. 


In  theory,  the  vectors  v,  computed  from  stage  1  of  this  process  form  an  orthonormal  basis  of 
the  Krylov  subspace  Km  —  span{ro,Aro,..Am~lro}.  It  can  be  easily  verified  that 

A0^m  —  VmTm  +  0m+lvm+l^m  (3.7) 

and  therefore  that  V^AVm  =  Tm  which  means  that  Tm  is  nothing  but  the  matrix  representation 
of  the  section  of  A  in  the  Krylov  subpace  Km  with  respect  to  the  basis  Vm.  Furthermore,  it  is 
easily  seen  that  the  Lanczos  algorithm  realizes  a  projection  process,  i.e.,  a  Galerkin  process,  onto 
the  Krylov  subspace  Km  [13,  16].  In  other  words  the  approximate  solution  can  be  found  by 
expressing  that  it  belongs  to  the  affine  subspace  +  Km  and  that  its  residual  vector  b  -  is 

orthogonal  to  Km.  Denoting  by  IIm  the  orthogonal  projector  onto  I(m,  this  means  that  the  Lanczos 
method  solves  the  approximate  problem: 

Find  x(m)  €  x(°)  +  Km  such  that  : 

nm(b-  Ax(ro>)  =  0  (3.8) 

The  approximation  thus  computed  is  identical  with  that  provided  by  m  steps  of  the  conjugate 
gradient  (CG)  method  when  A  is  positive  definite  [13].  When  A  is  not  positive  definite  this  relation 
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between  the  Lanczos  algorithm  and  the  CG  method  can  be  exploited  to  derive  stable  generalizations 
of  the  CG  algorithm  to  symmetric  indefinite  systems  [12,  13,  6,  17).  In  this  paper  we  will  show 
smother  way  of  exploiting  this  relation  to  provide  a  method  for  treating  nearly  singular  systems, 
which  can  be  viewed  as  a  variation  of  the  conjugate  gradient  method. 

4.  The  Lanczos  Deflation  Algorithm 

According  to  the  previous  section,  when  A  is  not  nearly  singular  the  approximate  solution  to 
(1.1)  can  be  obtained  by  a  Lanczos  •  conjugate  gradient  method  with  the  solution  given  by: 


x(*)  =  x«»+Vt2<*> 

where  z^  is  a  Ar-dimensional  vector  which  satisfies  the  equation 

Tk*W  =  ||ro||ei.  (4.1) 


When  A  is  nearly  singular,  a  difficulty  arises  in  the  solution  of  the  above  tridiagonal  system. 
Indeed  for  large  £  it  is  well-known  that  the  extreme  eigenvalues  of  Tk  will  converge  to  the  corre¬ 
sponding  extreme  eigenvalues  of  A,  see  e.g.  [14],  and  as  a  result  the  eigenvalues  closest  to  zero  of 
the  matrix  Tk  for  large  enough  k  will  be  close  to  the  eigenvalues  closest  to  zero  of  A,  i.e.,  it  will  be 
just  as  nearly  singular. 

The  obvious  remedy  is  to  solve  (4.1)  with  either  of  the  two  deflation  procedures  of  Section  2. 
Let  us  assume  that  we  compute  the  smallest  eigenvalue  \[k *  and  the  associated  eigenvector  of 
the  tridiagonal  matrix  Tk  by  inverse  iteration.  If  we  apply  the  deflation  techniques  of  Section  2  to 
the  tridiagonal  system  (4.1),  we  obtain  the  decomposition 


+  2 —Jk) 
+  \l*>  1 


(4.2) 


where  A^  is  the  smallest  eigenvalue  of  Tk,  u[k^  is  the  eigenvector  associated  with  A|*\  z^  is  the 
deflated  solution  of  (4.1)  and 

i(k)  =  lko||«[«it). 

Multiplying  both  members  of  the  above  equation  on  the  left  by  the  matrix  Vk  and  adding  *(°)  we 
obtain 

x(*)  =  z(°)  +  V**(*)  =  *(°)  +  Vkz{k)  +  ^Vku[k).  (4.3) 

* l 

Observe  that  in  the  above  equation  the  vector  Vku[k ^  is  nothing  but  the  Ritz  vector  associated 
with  the  eigenvalue  of  A  closest  to  zero,  i.e.,  it  is  the  approximate  eigenvector  of  A  computed  by 
the  Lanczos  process  from  the  Krylov  subspace  Iik.  We  will  denote  it  by  w[kK  Hence  we  obtain  am 
approximate  deflated  solution  in  the  form 


where 


(*) 


,(*)  _  ,(*>  + 

X  -xd  +  >(i)u/i  , 


x{k)  =  xf°)  +  vkz{k). 


(4.4) 


(4.5) 


We  point  out  the  important  fact  that  (4.4)  is  the  orthogonal  decomposition  of  the  vector 
XW  -  a;(°)  in  the  direction  w[k '  in  the  subspace  A'*.  This  is  true  because,  as  may  be  readily  shown, 
VieZj ^  is  orthogonal  to  the  Ritz  vector  wj**. 

Note  that  we  need  to  store  the  vectors  r,  as  they  are  generated  and  retrieve  them  once  when 
forming  the  approximation 

5.  Practicalities 

5.1.  Loss  of  Orthogonality  of  the  Lanczos  Vectors 

A  troublesome  behavior  of  the  Lanczos  algorithm  is  the  loss  of  orthogonality  of  the  vectors 
v i,t  =  l,...m.  Fortunately,  this  does  not  prevent  the  method  from  converging  but  often  results 
in  a  slow-down.  Parlett,  Scott  and  Simon  [13,  18,  19]  have  proposed  several  different  practical 
reorthogonalization  techniques.  Loss  of  orthogonality  is  a  phenomenon  that  cannot  be  avoided 
without  some  form  of  reorthogonalization  but  can  be  delayed  by  replacing  the  trivial  implementation 
(3.2)  -  (3.5)  by  the  following  one  [14,  11,  18] 

q  :=  Avj  - 
<*j  ■= 

q:=q-  ajvj. 

Then  Vj+\  =  q  and  (3.4)  and  (3.5)  deliver  the  next  vector  Vj+\.  This  simply  corresponds  to 
a  modified  Gram-Schimdt  step,  instead  of  the  regular  Gram-Schmidt  method  of  (3.2)  -  (3.4). 
Additional  reorthogonalization  can  be  added  to  further  post-pone  loss  of  orthogonality  at  little 
extra  cost,  by  appending  the  following  reorthogonalization  steps  to  the  above. 

6  :=  (g,  v/-i) 
q  ••=  q  -  6vJ_l 
6  :=  (g,  vj) 

otj  :=  atj  +  6 

g  :=  g  -  Svj 

Again  define  vy+i  :=  q  and  apply  (3.4)  and  (3.5)  to  get  v/+i. 

So  far  we  have  not  discussed  the  possible  negative  effects  of  the  loss  of  orthogonality  in  the 
Lanczos  process  mentioned  earlier.  An  important  factor  of  this  phenomenon  is  the  way  in  which 
it  appears.  Basically,  loss  of  orthogonality  is  a  signal  that  one  or  a  few  approximate  eigenvalues 
of  the  matrix  A  have  reappeared  after  they  have  already  converged.  As  a  result  one  can  expect 
that  as  soon  as  the  eigenvalue  A^  has  converged  to  Aj  then  a  second  copy  of  this  eigenvalue  will 
appear  in  the  following  steps  of  Lanczos.  The  presence  of  this  extra  eigenvalue  can  be  disastrous 
if  we  solve  the  tridiagonal  system  without  some  extra  precaution,  because  we  are  now  facing  a 
tridiagonal  system  with  multiple  sigularity. 

As  will  be  seen  in  Section  6  the  first  eigenvector  will  likely  converge  at  the  same  time  as  the 
deflated  solution  converges,  so  this  phenomenon  will  seldom  hamper  the  progress  of  our  procedure. 
However,  as  is  discussed  in  section  5.2  a  reasonable  computational  code  should  foresee  hard  cases 
that  might  occur  such  as  when  the  smallest  eigenvalue  of  A  is  so  well  separated  from  the  others 
that  convergence  of  the  eigenvalue  is  very  fast.  There  are  two  possible  remedies  to  the  problem 
of  loss  of  orthogonality.  The  first  one  is  to  simply  deflate  more  carefully  in  the  tridiagonal  system 
solution,  i.e.,  to  deflate  by  as  many  eigenvectors  as  there  are  small  eigenvalues. 
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The  second  solution,  which  is  well  knowm  in  the  context  of  the  eigenvalue  problem  [14,  18)  is 
to  perform  a  selective  reorthogonalization  (SO)  in  the  Lanczos  process.  Here,  in  contrast  with  full 
reorthogonalization,  one  only  reorthogonalizes  the  current  Lanczos  vectors  against  the  eigenvectors 
that  have  converged.  Moreover,  a  clever  implementation  allows  one  to  orthogonalize  only  when 
necessary.  The  idea  of  selective  reorthogonalization  is  intuitively  simple:  since  it  is  known  that  loss 
of  orthogonality  appears  mainly  in  the  direction  of  the  converged  Ritz  vectors,  then  a  remedy  is  to 
remove  the  corresponding  components  from  the  Lanczos  vectors  as  soon  as  these  start  to  reappear. 
The  important  point  is  that  there  are  ways  of  determining  when  these  components  are  likely  to 
come  back,  for  more  details  see  [14,  19). 

The  above  discussion  suggests  in  fact  that  we  may  only  have  to  reorthogonalize  v,-  against  the 
converged  Ritz  vector  that  corresponds  to  the  eigenvalue  closest  to  zero. 

5.2.  Stopping  Criteria 

When  implementing  the  algorithms  presented  in  this  paper  into  a  software  package,  one  must 
design  stopping  criterion  for  the  Lanczos  iteration.  There  are  two  independent  factors  affecting  the 
stopping  criterion.  First,  we  must  be  sure  that  the  eigenpair  (Aj*\  u^)  of  7*  has  already  converged 
to  the  eigenpair  (Aj,  inj)  of  A,  as  it  is  only  reasonable  to  compute  the  deflated  solution  of  Tk  after 
this  has  occurred.  Secondly,  we  should  stop  iterating  when  \\xd  -  x^  ||  is  reasonably  small. 

A  well-known  result  in  the  Lanczos  algorithm  is  that  it  is  not  necessary  to  compute  the  Ritz 
vector  in  order  to  check  for  convergence.  This  is  because  of  the  very  useful  relation  (14) 

AinJ*1  =  x\k)w\k)  +  pk+lelu\k\  (5.1) 

from  which  one  derives  the  residual  norm 

|«A  -  =  0k+l\elul%  (5.2) 

In  other  words  the  residual  norm  of  the  Ritz  pair  a)*\  is  equal  to  the  absolute  value  of  the  last 

component  of  the  eigenvector  uj^  of  the  tridiagonal  matrix  multiplied  by  0k+i .  This  provides  an 
inexpensive  way  of  checking  the  convergence  of  the  eigenpair  since  the  error  on  the  eigenvalue  is  of 
the  order  of  the  square  of  the  residual  norm  in  (5.2)  while  the  angle  between  the  exact  eigenvector 
and  the  approximate  one  is  of  the  same  order  as  the  residual  norm  [14]. 

A  formula  similar  to  (5.2)  can  be  easily  derived  for  the  residual  of  the  approximate  solution 
x^K  In  fact  we  must  first  define  what  might  be  a  suitable  analogue  to  the  residual  norm  in  the 
context  of  nearly  singular  systems.  Ideally,  we  would  like  to  consider  the  residual  norm  for  the 
matrix  PA  =  (I  -  wiwf)A,  i.e.,  we  would  wish  to  consider  P(b  -  Axj),  the  deflated  residual,  as  an 
appropriate  analogue  of  the  usual  residual  vector.  The  reason  why  this  is  the  correct  analogue  of  the 
residual  norm  in  the  non-sigular  case  is  that,  with  respect  to  xd,  we  are  in  fact  attempting  to  solve 
the  system  PAxd  =  Pb  in  the  subspace  PRn.  Unfortunately,  the  exact  projector  P  depends  on  the 
exact  eigenvector  u/j  which  is  not  known  a  priori.  However,  once  the  approximate  eigenvector 
has  converged,  we  can  use  it  in  place  of  wi  and  therefore  define  the  approximate  projected  residual 

r(dk)  =  pro (6  -  Ax{k))  =  (/- wj*1  [4M] T)  {b  ~  Ax[k)),  (5.3) 

where  pW  is  the  projector  in  the  direction  orthogonal  to  w[k\  Substituting  equation  (4.5)  in  the 
above  equation  we  get 

r<‘>  =  P<‘>  (b-A  [*<°>  +  V**<*>])  =  P<*>  (r<°>  -  AVkz^)  . 
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Using  the  relation  (3.7)  and  recalling  that  v\  =  ro/||ro||,  we  obtain 

rW  =  p(k)  ^||ro||Vl  _  AVkz{dk])  =  PW  (llmllvj  -  VkTkz{k)  -  /W^Wi) 

=  P{k)Vk  (ikolki  -  Tkz[k))  -  0k+leTkz{k)vk+ 1 

Now  observe  that  by  (4.1)  and  (4.2)  the  term  ||roj|ei  -  Tkz ^  is  equal  to  a  scalar  multiple  of  u^k\ 
and  since  Vku\k^  =  u>j**  the  first  term  in  the  above  sum  vanishes.  Hence  we  have  proved  that  the 
residual  norm  is  given  by 

A+i|e**?)N  (5-4) 

which  can  be  computed  at  little  cost  at  each  iteration. 

Peculiar  situations  may  arise  in  which  the  convergence  of  xd  ,  as  measured  by  the  above  resid¬ 
ual  norm,  occurs  before  the  eigenpair  of  Tk  has  converged.  Should  this  happen,  the  corresponding 
solution  x<*>  must  not  be  accepted.  Accordingly,  a  general  strategy  might  proceed  as  follows.  We 
iterate  with  the  Lanczos  process,  without  computing  z ^  or  x^d  \  until  the  estimate  (5.2)  indicates 

that  the  eigenpair  (A^,u;[*^)  has  converged.  The  eigenvector  that  is  used  in  (5.2)  can  be 
computed  every  few  iterations  by  inverse  iteration.  After  this  eigenpair  has  converged,  we  performs 
selective  orthogonalization  with  respect  to  this  Ritz  pair  in  order  to  prevent  it  from  reappearing 

/  JL\ 

in  Tk.  At  the  same  time,  we  start  computing  xd  either  explicitly  by  one  of  the  two  deflation 
algorithms  outlined  in  Section  2,  or  by  an  updating  procedure  similar  to  one  used  in  SYMMLQ 
[12].  When  the  estimate  (5.4)  indicates  that  xj,**  has  converged,  we  stop  the  Lanczos  iteration. 

6.  Theoretical  error  bounds 

In  this  section  we  address  the  issue  of  convergence  rates  and  will  derive  theoretical  error 

bounds  on  the  approximate  deflated  solution  x^ .  For  this  purpose  we  need  some  additional 

notation.  Recall  that  11*  is  the  orthogonal  projector  onto  the  Krylov  subspace  Kk.  We  denote  by 
qW  =  p(k) n*  the  orthognal  projector  onto  the  subspace  Iik  of  I(k ,  which  is  orthogonal  to  the 
Ritz  vector  w[k\  Thus,  the  rank  of  II*  is  k  in  general  while  the  rank  of  Q^k\  i.e.,  the  dimension  of 
the  subspace  Kk  =  QWRn  =  &k>*Kk,  is  k  -  1  in  general.  Without  loss  of  generality  it  is  assumed 
throughout  this  section  that  x^0)  =  0,  i.e.,  r^0)  =  b.  We  will  establish  our  result  with  the  help  of  a 
few  simple  lemmas. 

Lemma  6.1.  The  Ritz  deflated  solution  x^  is  the  (unique)  solution  of  the  Galerkin  problem 

Q{k)(b-  Ax)  =  0,  x  eQ{i) RN,  (6.1) 

Proof.  From  the  comments  following  (4.4),  (4.5),  it  is  clear  that  x^  belongs  to  the  subspace 

qWrW .  We  derive  from  the  Lanczos  relation  (3.7)  that 

Q{k)(b  -  Ax {k))  =  QW(b  -  AVkz(k))  =  Q<‘>  (||6||v,  -  V*r*4*>  -  0k+1vk+1el z[dk)) 

The  term  Q^vk+i  vanishes  because  vk+i  which  is  orthogonal  to  the  subspace  Kk,  is  also  orthogonal 
to  the  subspace  QWRn  of  Kk.  We  are  left  with 

Qlk)(b  -  Ax{k))  =  Q<*>  (g(t)||6||v,  -  VkTkz(k)) 
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Noticing  that  Q^v\  =  v\  -  vfw[k^w[k^  =  Vk(e\  -  efu^u^)  we  finally  obtain 
Q{k)(b  ~  Ax[k])  =  Q^Vk  (||6||(ei  -  efu^u^)  -  Tkz(k)) 

By  definition  of  the  deflated  solution  z\f\  the  above  vector  is  zero. 

I 

We  now  assume  that  A  is  positive  definite  and  denote  by  ||.||a  the  A-norm,  i.e.,  ||x||A  = 
(Ax,  x)1/2.  In  classical  Galerkin  techniques,  the  right  hand  side  b  in  (6.1)  belongs  to  the  subspace 
of  projection  so  that  we  approximately  solve  Ax  =  b  in  that  subspace.  Here,  however,  the  right 
hand  side  is  not  in  the  subspace  of  projection.  Observing  that  (6.1)  also  reads  {Q^ b  —  Ax)  =  0, 

we  can  say  that  the  Lanczos  deflation  method  attempts  to  solve  the  linear  system 

Ax  =  Q(*>6  (6.2) 

by  a  Galerkin  process  onto  the  subspace  Kk.  Thus,  for  the  approximation  to  be  accurate,  we  must 
show  that  the  projected  right  hand  side  Q^b  is  in  some  sense  close  to  the  deflated  right  hand 
side  Pb.  This  will  be  considered  in  detail  later.  We  now  apply  a  classical  argument  in  Galerkin 
methods. 

In  the  following  discussion  we  denote  by  xk  the  exact  solution  of  the  linear  system  (6.2).  As 
stated  above,  x^  is  the  Galerkin  solution  of  the  new  linear  system  (6.2),  i.e.,  it  is  a  Galerkin 
approximation  to  x*  from  the  subpsace  Kk.  A  standard  theorem  relating  the  Galerkin  method  to 
the  Rayleigh-  Ritz  method  yields  the  following  result  on  the  error  x^  -  x*. 

Lemma  6.2.  The  approximate  deflated  solution  x^  minimizes  the  function 

J{x)  =  ||x-  xk\\A 

among  all  elements  x  of  the  subspace  Kk  =  QWRN . 

We  now  wish  to  reformulate  the  above  results  in  terms  of  polynomials.  Clearly,  a  vector  v  is 
in  the  Krylov  subspace  Kk  if  and  only  if  it  can  be  expressed  as  v  =  p(A)b  where  p  is  a  polynomial 
of  degree  <  k  —  1.  The  next  lemma  shows  how  the  additional  constraint  that  v  belongs  to 
translates  for  the  polynomial  p. 

Lemma  6.3.  A  vector  v  oPRN  belongs  to  RN  if  and  only  if  it  is  of  the  form  v  =  q[A)b  where  q 
is  a  polynomial  of  degree  <  k  -  1  such  that  q(A\  ')  =  0. 

Proof.  Consider  a  vector  of  Kk  of  the  form  v  =  q(A)v\,  where  q  is  of  degree  <  k  —  1.  This  vector 
is  in  Kk  if  and  only  if  (v,wj**)  =  0,  i.e.,  if  and  only  if 

(q(A)vi,w[k))  =  0.  (6.3) 

The  section  of  the  linear  operator  A  in  A'*,  i.e.,  the  rank  k  linear  operator  A*  =  TlkA\Kk  ap¬ 
proximates  A  in  the  subspace  Kk  in  the  Galerkin  process:  its  matrix  representation  in  the  basis 
V*  is  VkTkVk  .  It  is  easy  to  show  that  AJkv\  =  UkA}vi  for  any  j,j  <  k  and  therefore  we  have 
q[A)v i  =  II*g(A*)vi,  which  substituted  in  (6.3)  yields 


0  =  (Tltq(At)vuw[k))  =  (vi,q(Ak)Ukw[k))  =  (v1,?(A*)u;j<:)). 


(6.4) 


Clearly,  w[k ^  is  an  eigenvector  of  At  associated  with  the  eigenvalue  and  as  a  result  the  above 
relation  becomes 

(vi,g(A(1fc))w[t))  =  q(x[k)){vuw[k})  =  0  . 

Notice  that  the  inner  product  (t/i,u;ifc')  is  the  first  component  of  the  eigenvector  of  the  tridiag- 

(Jb) 

onal  matrix,  associated  with  the  eigenvalue  X\  1 .  This  component  cannot  be  zero  for  a  nonreducible 
tridiagonal  matrix,  so  we  conclude  that  v  belongs  to  if*  if  and  only  if  =  0  . 

I 

Denoting  by  Pk  the  space  of  all  polynomials  q  of  degree  <  k ,  such  that  g(A^)  =  0,  we  can 
state  an  immediate  corollary  of  the  above  lemmas. 

Corollary  6.1.  The  approximation  is  such  that 

\\xdk)  -  ZfclU  =  min  h(A)b  -  A_1Q(fc)6|U. 


In  the  next  proposition  we  use  this  equality  to  estimate  the  distance  between  the  exact  solution 
it  of  (6.2)  and  its  Galerkin  approximation  x^. 

Proposition  6.1.  Assume  that  k  is  sufficiently  large  that  —  A^  >  Ai.  Then  at  the  k-th  step  of 
the  deflated  Lanczos  procedure  we  have  the  inequality 


I \*dk)  ~  **IU  5  AlA(*)Alc'tHbllA  +  (^7  +  ««e(wi»wifc))IWI  + 

where  Tk-i  is  the  Chebyshev  polynomial  of  degree  k  -  1  of  the  Brst  kind. 


l  +  k*{l  +  vk) 


(6.5) 


vk  = 


Ay  +  Aa- A[fc) 
Aw  -  Aa  +  A|4)  ’ 


Ct/k  converges  to  ,n  which  v  =  lim  vk,  and  0(r,  y)  represents  the  acute  angle  between 

the  vectors  x  and  y  in  RN. 

Proof.  Recall  that  P  denotes  the  (eigen)-projector  onto  the  subspace  orthogonal  to  uq,  i.e.,  P  = 
I  -  wivt/i-  For  any  polynomial  q  in  Pk  we  have 


||*ifc)  -  xk\\A  <  lk(A)6  -  A-'QWb\\A  <  ||?(A)6  -  A~lPb\\A  +  llA^P  -  Q<‘>)t|U 

<  \\q(A)b  -  A~lPb\\A  +  ||A“1P(P  -  $<*>)6|U  +  ||A-‘(/  -  P)(P  -  QW)6||a 


or, 


||x'fc)  -  xk\\A  <  |k(A)6  -  A~lPb\U  +  ||A-1P(/  -  Q(k))b\\A  +  HA"1  (I  -  P)Q^b\\A.  (6.6) 

Consider  first  the  last  two  terms  of  the  right  hand  side.  Since  Q^b  =  p(k>lb,  by  using  elementary 
properties  of  projectors,  we  get 


||A_lP(J  -  Q^)b\\A  =  || A~lPP(I  -  Pw)b\\A  <  ||A-1P||a||P(J  -  P{k))b\\A 
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It  is  easy  to  show  that 

||A-‘P|U  =  ||R4-‘|U  =  1  (6.7) 

Moreover,  ___ 

||P(/-  p(‘i)6|U  <  V^||P(/-  pW)b\\  <  v^||P(/~  p(t))||||6(|, 

and 

||P(/-P(*>)|l  =  5.»e(tt;1,«;lfc)).  (6.8) 

The  last  term  in  the  right  hand  side  of  (6.6)  satisfies 

IIA-'V-  P)g^6|U  =  II A-1  (I  -  P)(I  -  ^P^IU  <  IIA-^J-  P)|U||(/  -  P)P(*)6|U, 

with  ||A-1(/-  P)||A  =  1/Ai.  The  vector  (I  -  P)Q^b  ,  when  expanded  with  respect  to  the 
eigenvectors  of  A,  has  only  one  term  corresponding  to  the  first  eigenvector.  Hence, 

IIU  -  P)P{k]b\u  =  \ATllU  -  P)pWb\\  <  n/X7||(j  -  p)p<*i||||b||, 

and  it  is  again  possible  to  show  that 

||(/_  P)p(*)|j  =  Sine(wi,w[k)).  (6.9) 

Using  these  upper  bounds  for  the  last  two  terms  of  (6.6)  we  arrive  at 

II -  *fclU  <  II q(A)b  -  A_lP6|U  +  stn©(u>1,u'|i!))||6||.  (6.10) 

We  now  seek  a  particular  polynomial  q  in  Pk  for  which  the  first  term  in  the  above  expression 
is  as  small  as  possible.  Consider  the  particular  polynomial  of  degree  k  defined  by 

p(A) = 7w  [Af(A  -  Ait))  -  (a  -  . 

where  f(A)  is  the  polynomial  of  degree  k  —  1  defined  by 

,/u  _  Tk-i{vk  ~  at  A) 

in  which 

__  Ajy  +  A2  —  _  1  +  i/k  __  2 

"““a,,- *,  +  ><*»’  °‘~~~A,v-Aj  +  a!  *»' 

The  polynomial  t(A)  is  of  degree  k  -  1  and  its  value  at  the  origin  is  1.  In  fact,  it  is  chosen  so  as 
to  minimize  the  infinity  norm  in  the  interval  [A2  -  a[*\  A/v],  over  all  such  polynomials.  We  observe 
that  p(Aj*))  =  p(0)  =  1  which  means  that  p  can  be  written  as  p(A)  =  1  —  Ag(A).  Moreover,  q  is  of 
degree  k  -  1  and  we  have  ?(Aj^)  =  0  so  that  the  polynomial  q  is  in  I(k. 

Let  us  expand  the  vector  A-1f>  in  the  eigenbasis  of  A  as  A-16  =  YLiLi  7iu,i-  Noticing  that 
PA~lb  =  A~lb  -  'iiwi  we  get 

lk(A)6  -  A-lP6|ft  =  ||g(A)AA-‘6  -  PA~lb\\\  =  A,[l  -  p(A1)]27?  +  £  A  ,-7?p(A,)2.  (6.11) 

1-2 


We  start  by  focussing  on  the  first  term  of  the  right  hand  side  in  the  above  equality.  The  factor 
p(Ai)  —  1  satisfies  the  relations: 

p(a.)  - 1  =  -jSfifAi  -  a',*’)  -  ^-<s>  l(Al)  -  : 

Aj 

-  4w»i  -  M*’>  - 1)  +  - 1)- 

At  At 

Remembering  that  <(0)  =  1  this  can  be  rewritten  as 

p(Aj)  -  1  Ax  -  A<*>  f <(Ax  -  \[k))  -  «(0)1  ,  A(.fc)  -  Ax  r<(Ax)  -  t(0) 


i  Ax 

where  Ax  —  A^  <  fx  <  0,  0  <  &  <  Ax.  The  derivative  t'{ A)  is  given  by 


t'( A)  =  -(k  -  l)a* 


Uk-ijvk  ~  at  A) 

Tk-\{vk) 


(6.12) 


where  Uk-2  is  the  Chebyshev  polynomial  of  degree  k  -  2  of  the  second  kind.  Since  Ax  -  a|** 
converges  to  zero  from  the  left,  and  since  the  function  \Uk-2{vk  -  a*A)|  decreases  in  the  interval 

(Ax  -  Ai*\Ax],  we  have  |</(f2)|  <  (^(Ci)!  <  (t'(Ax  -  a(^)|  s  Ck  .  It  is  clear  that  at  the  limit  <7*  is 
equivalent  to 

which  in  turn  can  be  shown  to  be  asymptotic  with 


/•T-1 


Finally,  going  back  to  the  first  term  in  (6.11), 


A.[1  -p(a,)I!7?  =  <  iL-diilpmuj  < 

*  *  A  j 

Using  the  inequality  (o2  +  62)1/2  <  |a|  +  |6|  in  (6.11)  and  the  above  bound  we  obtain 

A(*)_Al  11/2 
\\q(A)b  -  A-'Pb\\A  <  2Ck  1  (<t)  l\\b\\A  +  £  AnfaA,)2 

*1  Li-2 

Consider  now  the  second  part  of  the  right-hand-side  in  (6.13).  We  have 

P(A.)  =  A.-  [‘(AiL-  <(A;  -  A}*1)]  +  <(A.  .  Aj*lj 


(6.13) 


(6.14) 
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By  the  mean-value  theorem,  the  expression  between  brackets  can  be  expressed  as  £'(£),  the  deriva¬ 
tive  of  t(A)  at  some  point  £  in  the  interval  [A,-  -  a(*\a,].  The  point  £  belongs  to  the  interval 

[A2  —  Ai*\Ajy]  and  therefore  its  transformed  value  {vk  -  otkO/^k  belongs  to  the  interval  [-1,1]. 
From  the  expression  (6.12)  and  the  fact  that  |C7’*,_2 (x) |  <  (k  -  2)  for  x  €  [-1, 1],  it  is  clear  that 
t'{€)  <  k3ak/Tk- i(^i).  Moreover,  the  second  term  of  of  the  right-hand  side  of  (6.14)  is  naturally 
bounded  by  l/Tk-i{uk).  Therefore, 


Ip(A,)|  <  A.; 


k* 


+ 


1 _  1  +  k2ak\N  _  1  +  &2(1  +  vk) 

Tk-i(vk)  ~  Tk-i(vk)  Tk-iivk) 

Thus,  the  expression  between  brackets  in  (6.13)  is  bounded  from  above  as  follows 


N 


E  Vyfr(*»)s 


L«=2 


1/2 


1  +  fc2(l  +  ISk) 

Tk-iWk) 


N 


EW 


Li=2 


1  +  k2(l  +  vk) 

Tk-xWk) 


IMU- 


(6.15) 


(6.16) 


The  result  follows  by  combining  (6.16),  (6.13),  and  (6.10). 


I 


A  few  comments  on  the  above  proposition  are  in  order.  The  term  sin  Q(wi,w 
to  zero  like  [14] 


!fc)) 


1 

TthiV 


converges 

(6.17) 


where 


'll 


Ay  +  A2  —  At 
A n  -  A2 


(6.18) 


Similarly,  the  relative  error  on  the. eigenvalue  (A^  -  Aj)/A^  converge  to  zero  decreasingly  as  the 
square  of  the  quantity  (6.17).  Observe  the  similarity  between  the  number  71  and  the  coefficient  uk 
of  the  proposition.  Since  A^  is  close  to  zero  for  large  k,  these  two  numbers  will  be  close  to  each 
other  at  the  limit.  The  coefficient  Ck  in  the  proposition  is  not  bounded  but  is  of  the  order  of  O(k). 

However,  its  product  with  (Aj^  -  Ai)/A^  tends  to  zero  rapidly.  The  same  observation  holds  for 
the  term  k3  in  the  numerator  of  the  last  term  of  (6.5). 

We  now  turn  our  attention  to  the  actual  error  e*  =  xj  —  xjfK  We  cannot  use  an  A— norm  to 
measure  this  error  because  if  Aj  is  small  this  norm  will  dampen  the  component  of  the  error  in  the 
wi  direction.  As  a  result  a  possible  large  component  in  that  direction  can  be  unfairly  hidden  by 
this  norm.  Therefore,  we  choose  to  use  the  usual  Euclidean  norm  ||.||.  Moreover,  we  separate  the 
above  error  in  two  distinct  parts,  namely  the  component  (/  -  P)ek  in  the  uq  direction,  and  the 
component  Pek  orthogonal  to  it.  We  show  that  both  terms  tend  to  zero  quickly  as  Jfc  increases. 


Proposition  6.2.  The  part  of  the  error  in  the  direction  of  the  eigenvector  w\  satisfies  the  inequality: 


<  sin  ) 


(6.19) 


Proof.  We  have 

II U  -  P)(*.  -  4*')ll  -  IIU  -  r)4*’«  =  -  J’)PW4‘)II  <  lltf  -  P)P“,llll4‘,ll 


The  result  follows  from  (6.9). 

I 

The  above  result  means  that  relatively  to  ||x^||,  the  error  in  the  direction  of  w\  is  bounded  by 
the  sine  of  the  angle  between  the  exact  eigenvector  and  the  Ritz  vector.  As  we  mentioned  above, 
this  angle  is  known  to  converge  to  zero  as  rapidly  as  the  sequence  (6.17)  see  [14]. 
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Proposition  6.3.  The  part  of  the  error  in  the  eigeuspace  orthogonal  to  wi  satisfies  the  inequality: 

\\P[xd  -  4*,)||  <  e(u>i,u>[*,)||6||  +  ~^=,  (6.20) 

where  e*  is  the  right  hand  side  of  the  inequality  (6.5). 

Proof.  By  the  triangle  inequality, 

ii  p{*4  -  xifc))ii  <  ii  p(xd  -  it)n  +  up(4fc)  -  £*)ii. 

Since  A~lP  =  PA~l  =  PA~lP ,  the  first  term  of  the  right  hand  side  is  such  that 
|| P(xd  -  it) ||  ==  I! P(A~lPb  -  A~l  P<*>6)||  =  ||PA-‘(P6  —  P<*>6)|| 

=  |!A"lPP(/-P<*))6||  <  HA^PjlllPCf  -  P(t)6)||||6||  <  -Utn  ©(wi,  w}*>)||6||. 

For  the  second  term,  using  the  inequality 

||Py||  <  — |=||Py|U, 

we  get  immediatly  that 

intf1  -  «)ii  <  -  »>iu  <  -  *tiix  s 

I 

7.  Numerical  Experiments 

We  now  present  the  results  of  some  numerical  experiments  to  verify  the  accuracy  and  stability 
of  the  Lanczos-deflation  algorithms.  All  computations  were  performed  on  a  VAX-780  in  single 
precision,  with  a  relative  machine  precision  of  about  10~7. 

In  these  experiments,  we  employ  a  simple  stopping  criterion.  The  convergence  of  the  eigenpair 
(A^, Uj*')  is  checked  every  five  Lanczos  iterations.  The  tolerance  in  the  stopping  criterion  is  set 
to  10~5.  When  this  pair  has  converged,  we  stop  the  algorithm.  No  further  test  on  the  convergence 
of  xjp  is  made.  In  fact,  in  all  of  our  tests  the  convergence  of  x^  occurs  almost  simultaneously. 

In  a  first  test  we  solve  the  linear  system  AiX  =  b  where  Ai  =  d*ap(10_/,2,3,  •  •  •  ,n),  with  I 
varying  from  1  to  8;  bT  =  (1,  •  •  • ,  1),  and  n  =  100.  The  deflated  solution  to  the  above  problem 
is  xd  =  (0, 1/2,  ■•,  1/i,  ••,  1/n).  We  computed  the  deflated  solution  using  the  Lanczos-Projection, 
and  Lanczos-QL  methods.  For  comparison,  we  also  used  a  standard  conjugate  gradient  method  to 
solve  the  test  problem  and  then  deflated  the  solution  as  follows  : 

xd  =  x  -  (xTui)ui,  where  uf  =  (1,0,  •  •  •  ,0). 

Figure  1  shows  the  relative  error  in  xd  versus  I  for  the  three  methods.  Figure  3  illustrates  the 
simultaneous  convergence  of  the  approximate  eigenpair  and  deflated  solutu  i. 

Our  second  test  repeats  the  above  experiments  with  the  matrix  A<i  =  T  —  (Ai  —  <r )/,  where 
T  =  Tridiag{-1,2,  -1},  Ax  is  the  smallest  eigenvalue  of  T  and  <r  =  10-/,  with  I  again  varying  from 
1  to  8  and  n  =  20.  Note  that  the  smallest  eigenvalue  of  At  is  a  when  a  is  small.  The  solution  is 


chosen  such  that  x  =  x<j  +  wi  where  xj  is  such  (xj,  w\)  =0.  The  right  hand  side  6  is  then  obtained 
by  forming  Axj  +  Aw\.  Figure  2  shows  the  relative  error  in  zj  versus  /  for  the  three  methods.  In 
Figure  3  we  illustrate  the  simultaneous  convergence  of  the  eigenpair  and  of  the  deflated  solution 
for  the  first  test  matrix  Ai  with  a  =  10-4.  The  plot  shows  the  residual  norms  of  the  eigenpair  and 
of  the  deflated  solution  as  given  by  (5.2)  and  (5.4)  respectively.  Figure  4  is  a  similar  illustration  of 
this  simultaneous  convergence  for  the  matrix  Ai. 

The  numerical  results  show  that  both  Lanczos  deflation  methods  compute  the  deflated  solution 
with  accuracy  close  to  machine  precision  independent  of  the  singularity  of  the  matrix  A.  On  the 
other  hand,  the  conjugate  gradient  method  without  deflation  becomes  unstable  as  A  becomes  more 
singular,  especially  in  the  second  example.  The  accuracies  displayed  by  the  two  deflation  methods 
are  very  similar. 

8.  Concluding  Remarks 

We  have  presented  two  different  ways  of  extending  the  Lanczos  algorithm  to  solving  nearly 
singular  systems.  Both  methods  retain  the  advantages  of  the  classical  Lanczos-Conjugate  gradient 
procedure  in  that  they  access  the  matrix  A  only  in  the  form  of  matrix  by  vector  products.  The 
overhead  of  the  algorithm  over  regular  Lanczos  is  limited  to  deflating  a  tridiagonal  matrix  and 
is  negligible  compared  to  the  cost  of  the  overall  computation.  There  doesn’t  seem  to  be  much 
difference  in  the  performance  of  the  two  deflation  techniques  for  tridiagonal  matrices.  While  the 
QL  method  is  slightly  more  complicated  than  the  orthogonal  projection  method,  it  should  be  more 
robust  and  more  easily  extensible  to  higher  dimensional  nullity  problems. 

Although  we  have  presented  an  algorithm  which  requires  the  storage  of  the  Lanczos  vectors,  we 
should  emphasize  that  an  updating  version  similar  to  SYMMLQ  could  easily  be  derived.  Moreover, 
the  techniques  developed  here  can  in  principle  be  extended  to  handle  higher  dimensional  null  spaces 
and  nonsymmetric  problems. 


Acknowledgement:  The  authors  thank  Mr.  Herbert  Beke  for  his  assistance  throughout  this 
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Figure  2:  Relative  Error  versus  degree  of  singularity  of  three 
methods:  Regular  CG  (CGM),  Lauczos-Projection  (L-P)  and 
Lanczos=QL  (L-QL). 
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